using System;
using System.Collections;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.IO;
using igp.net.DataModel;

namespace igp.net.DataModel
{
    /// <summary>
    /// Represents a polynomial.
    /// </summary>
    public class igpPolynomial : IgpFunctionBase
    {
        /// <summary>
        /// Default constructor.
        /// </summary>
        public igpPolynomial()
        {
            Monomials = new ArrayList();
            DomainLower = -5;
            DomainHigher = 5;
            Limited = false;

            ItemName = "Polynomial ";
            TypeID = IgpItemType.Polynomial;
            LineWidth = 1;
            LineColor = Color.Green;
            TrendlineName = "";
        }

        /// <summary>
        /// Gets the list of contained monomials.
        /// </summary>
        public ArrayList Monomials { get; private set; }

        /// <summary>
        /// Gets/Sets a PointSet that will be used to calculate parameters of a specific trendline.
        /// </summary>
        public igpPointSet TrendlinePoints { get; set; }

        /// <summary>
        /// Specifies the interpolation mode.
        /// </summary>
        public IgpInterpolationMode TrendlineMode { get; set; }

        internal String TrendlineName { get; set; }

        /// <summary>
        /// Adds a new monomial.
        /// </summary>
        /// <param name="m">Multiplier</param>
        /// <param name="w">Exponent</param>
        /// <returns>Added monomial</returns>
        public IgpMonomial AddMonomial(double m, double w)
        {
            var temp = new IgpMonomial(m, w);
            Monomials.Add(temp);
            return temp;
        }

        /// <summary>
        /// Calculates values of this function for specified arguments.
        /// </summary>
        /// <param name="x">Argument</param>
        /// <returns>f(x)</returns>
        public override double Value(double x)
        {
            double result = 0;

            foreach (IgpMonomial mon in Monomials)
            {
                if (!mon.Limited)
                {
                    result += mon.Multiplier*Math.Pow(x, mon.Exponent);
                }
                else
                {
                    if (x >= mon.LowerXLimit && x < mon.HigherXLimit)
                    {
                        result += mon.Multiplier*Math.Pow(x, mon.Exponent);
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Calculates values of the derivative of this function for specified arguments.
        /// </summary>
        /// <param name="x">Argument</param>
        /// <returns>f'(x)</returns>
        public override double Derivative(double x)
        {
            double result = 0;

            foreach (IgpMonomial mon in Monomials)
            {
                if (!mon.Limited)
                {
                    if (mon.Exponent > 0)
                    {
                        result += mon.Exponent*mon.Multiplier*Math.Pow(x, mon.Exponent - 1);
                    }
                }
                else
                {
                    if (x >= mon.LowerXLimit && x < mon.HigherXLimit && mon.Exponent > 0)
                    {
                        result += mon.Exponent*mon.Multiplier*Math.Pow(x, mon.Exponent - 1);
                    }
                }
            }


            return result;
        }

        /// <summary>
        /// Calculates values of the indefinite integral of this function for specified arguments.
        /// Constant factor is assumed to be 0.
        /// </summary>
        /// <param name="x">Argument</param>
        /// <returns>Indefinite Integral(x)</returns>
        public override double Integral(double x)
        {
            double result = 0;

            foreach (IgpMonomial mon in Monomials)
            {
                if (!mon.Limited)
                {
                    result += mon.Multiplier/(mon.Exponent + 1)*Math.Pow(x, mon.Exponent + 1);
                }
                else
                {
                    if (x >= mon.LowerXLimit && x < mon.HigherXLimit)
                    {
                        result += mon.Multiplier/(mon.Exponent + 1)*Math.Pow(x, mon.Exponent + 1);
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Saves data to a TextWriter.
        /// </summary>
        /// <param name="sw">data output</param>
        public override void SaveData(TextWriter sw)
        {
            sw.WriteLine("<polynomial>");
            sw.WriteLine("ItemName = {0}", ItemName);
            sw.WriteLine("LineColor = {0}", LineColor.ToArgb());
            sw.WriteLine("LineWidth = {0}", LineWidth);
            sw.WriteLine("LineType = {0}", (int) LineType);

            sw.WriteLine("DomainMinimum = {0}", DomainLower);
            sw.WriteLine("DomainMaximum = {0}", DomainHigher);

            if (TrendlinePoints != null)
            {
                sw.WriteLine("TrendlineMode = {0}", (int) TrendlineMode);
                sw.WriteLine("TrendinePoints = {0}", TrendlinePoints.ItemName);
            }
            else
            {
                IEnumerator item = Monomials.GetEnumerator();

                item.Reset();
                while (item.MoveNext())
                {
                    var p = (IgpMonomial) item.Current;
                    sw.WriteLine("AddMonomial = {0}|{1}|{2}|{3}|{4}", p.Multiplier, p.Exponent, p.Limited, p.LowerXLimit,
                                 p.HigherXLimit);
                }
            }

            sw.WriteLine("IsVisible = {0}", IsVisible);
            sw.WriteLine("OnLegend = {0}", VisibleOnLegend);

            sw.WriteLine("</polynomial>");
        }

        /// <summary>
        /// Loads data from a TextReader.
        /// </summary>
        /// <param name="sr">data input</param>
        public override void ReadData(TextReader sr)
        {
            String Line;
            while ((Line = sr.ReadLine()) != "</polynomial>")
            {
                char[] sep = {'=', '|'};
                String[] cmd = Line.Split(sep, 6);

                if (cmd[0].StartsWith("ItemName"))
                {
                    ItemName = cmd[1].Substring(1, cmd[1].Length - 1);
                }

                if (cmd[0].StartsWith("LineColor"))
                {
                    LineColor = Color.FromArgb(int.Parse(cmd[1]));
                }
                if (cmd[0].StartsWith("LineWidth"))
                {
                    LineWidth = int.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("LineType"))
                {
                    LineType = (DashStyle) int.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("LineWidth"))
                {
                    LineWidth = int.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("DomainMinimum"))
                {
                    DomainLower = double.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("DomainMaximum"))
                {
                    DomainHigher = double.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("TrendlineMode"))
                {
                    TrendlineMode = (IgpInterpolationMode) int.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("TrendinePoints"))
                {
                    TrendlineName = cmd[1].Substring(1, cmd[1].Length - 1);
                }
                if (cmd[0].StartsWith("AddMonomial"))
                {
                    IgpMonomial mon = AddMonomial(double.Parse(cmd[1]), double.Parse(cmd[2]));
                    mon.Limited = bool.Parse(cmd[3]);
                    mon.LowerXLimit = double.Parse(cmd[4]);
                    mon.HigherXLimit = double.Parse(cmd[5]);
                }


                if (cmd[0].StartsWith("Visible"))
                {
                    IsVisible = bool.Parse(cmd[1]);
                }
                if (cmd[0].StartsWith("OnLegend"))
                {
                    VisibleOnLegend = bool.Parse(cmd[1]);
                }
            }
        }

        private void TrendlineLagrange(igpPointSet ptset)
        {
            Monomials.Clear();
            int ilosc = ptset.Points.Count;

            if (ilosc <= 1) return;

            var main = new Matrix.Matrix(ilosc, ilosc + 1);
            int i = 0;

            IEnumerator item = ptset.Points.GetEnumerator();
            item.Reset();
            item.MoveNext();

            double xm = ((IgpPoint) item.Current).X;

            double tminx = xm;
            double tmaxx = xm;

            item.Reset();

            while (item.MoveNext())
            {
                for (int p = 0; p < main.Cols - 1; p++)
                {
                    if (((IgpPoint) item.Current).X != 0 || main.Rows - p != 1)
                    {
                        main.Data[i, p] = Math.Pow(((IgpPoint) item.Current).X, main.Rows - p - 1);
                    }
                }

                main.Data[i, main.Cols - 1] = ((IgpPoint) item.Current).Y;
                i++;

                xm = ((IgpPoint) item.Current).X;
                if (tminx > xm) tminx = xm;
                if (tmaxx < xm) tmaxx = xm;
            }

            DomainLower = tminx;
            DomainHigher = tmaxx;

            Matrix.Matrix data = main.Gauss();

            for (int x = 0; x < ilosc; x++)
            {
                AddMonomial(data.Data[x, 0], main.Rows - x - 1);
            }
        }

        private void TrendlineLinear(igpPointSet ptset)
        {
            Monomials.Clear();

            double srx = 0, sry = 0, D = 0;
            double a = 0;

            IEnumerator item = ptset.Points.GetEnumerator();
            item.Reset();
            int k = 0;

            while (item.MoveNext())
            {
                srx += ((IgpPoint) item.Current).X;
                sry += ((IgpPoint) item.Current).Y;
                k++;
            }

            if (k > 0)
            {
                srx = srx/k;
                sry = sry/k;

                item.Reset();

                while (item.MoveNext())
                {
                    a += ((IgpPoint) item.Current).Y*(((IgpPoint) item.Current).X - srx);
                    D += Math.Pow(((IgpPoint) item.Current).X - srx, 2);
                }

                a = a/D;
                double b = sry - a*srx;

                AddMonomial(b, 0);
                AddMonomial(a, 1);
            }
        }

        private void TrendlineBSplineNormal(igpPointSet ptset)
        {
            Monomials.Clear();

            int ilosc = ptset.Points.Count;

            if (ilosc <= 2) return;

            IEnumerator item = ptset.Points.GetEnumerator();
            item.Reset();

            item.MoveNext();
            double xm = ((IgpPoint) item.Current).X;
            double tminx = xm;
            double tmaxx = xm;

            while (item.MoveNext())
            {
                xm = ((IgpPoint) item.Current).X;
                if (tminx > xm) tminx = xm;
                if (tmaxx < xm) tmaxx = xm;
            }

            DomainLower = tminx;
            DomainHigher = tmaxx;


            var main = new Matrix.Matrix(4*(ilosc - 1), 4*(ilosc - 1) + 1);
            Matrix.Matrix dane;

//iteracja warunku zgodnosci wartosci //////////////////////////////////////////

//przegladanie po punktach

            IgpPoint cp, np;

            int firstColumn = 0;
            int row = 0;

            item.Reset();

            item.MoveNext();

            cp = ((IgpPoint) item.Current);
            item.MoveNext();
            np = ((IgpPoint) item.Current);

            while (cp != null && np != null)
            {
                main.Data[row, firstColumn] = Math.Pow(cp.X, 3);
                main.Data[row, firstColumn + 1] = Math.Pow(cp.X, 2);
                main.Data[row, firstColumn + 2] = cp.X;
                main.Data[row, firstColumn + 3] = 1;

                main.Data[row + 1, firstColumn] = Math.Pow(np.X, 3);
                main.Data[row + 1, firstColumn + 1] = Math.Pow(np.X, 2);
                main.Data[row + 1, firstColumn + 2] = np.X;
                main.Data[row + 1, firstColumn + 3] = 1;

                //wartosci w punktach
                main.Data[row, 4*(ilosc - 1)] = cp.Y;
                main.Data[row + 1, 4*(ilosc - 1)] = np.Y;

                row += 2;
                firstColumn += 4;

                cp = ((IgpPoint) item.Current);

                if (item.MoveNext())
                {
                    np = ((IgpPoint) item.Current);
                }
                else
                {
                    np = null;
                }
            }

//iteracja warunku zgodnosci I i II pochodnej //////////////////////////////////

//przegladanie po punktach

            item.Reset();

            item.MoveNext();
            item.MoveNext();
            cp = ((IgpPoint) item.Current);
            item.MoveNext();
            np = ((IgpPoint) item.Current);

            firstColumn = 0;

            while (np != null)
            {
                // I pochodna
                main.Data[row, firstColumn] = 3*Math.Pow(cp.X, 2);
                main.Data[row, firstColumn + 1] = 2*cp.X;
                main.Data[row, firstColumn + 2] = 1;

                main.Data[row, firstColumn + 4] = -3*Math.Pow(cp.X, 2);
                main.Data[row, firstColumn + 5] = -2*cp.X;
                main.Data[row, firstColumn + 6] = -1;

                // II pochodna

                main.Data[row + 1, firstColumn] = 6*cp.X;
                main.Data[row + 1, firstColumn + 1] = 2;

                main.Data[row + 1, firstColumn + 4] = -6*cp.X;
                main.Data[row + 1, firstColumn + 5] = -2;

                row += 2;
                firstColumn += 4;

                cp = ((IgpPoint) item.Current);
                if (item.MoveNext())
                {
                    np = ((IgpPoint) item.Current);
                }
                else
                {
                    np = null;
                }
            }

// 2 rownania na spline naturalny

            main.Data[4*(ilosc - 1) - 2, 0] = 3*Math.Pow(tminx, 2);
            main.Data[4*(ilosc - 1) - 2, 1] = 2*tminx;
            main.Data[4*(ilosc - 1) - 2, 2] = 1;

            main.Data[4*(ilosc - 1) - 1, 4*(ilosc - 1) - 4] = 3*Math.Pow(tmaxx, 2);
            main.Data[4*(ilosc - 1) - 1, 4*(ilosc - 1) - 3] = 2*tmaxx;
            main.Data[4*(ilosc - 1) - 1, 4*(ilosc - 1) - 2] = 1;

            dane = main.Gauss();

            int lastread = 0;
            item.Reset();
            item.MoveNext();

            cp = ((IgpPoint) item.Current);
            item.MoveNext();
            np = ((IgpPoint) item.Current);

            while (np != null)
            {
                for (int y = 0; y < 4; y++)
                {
                    IgpMonomial mon = AddMonomial(dane.Data[lastread, 0], 3 - ((lastread)%4));
                    mon.LowerXLimit = Math.Min(cp.X, np.X);
                    mon.HigherXLimit = Math.Max(cp.X, np.X);
                    mon.Limited = true;
                    lastread++;
                }

                cp = ((IgpPoint) item.Current);
                if (item.MoveNext())
                {
                    np = ((IgpPoint) item.Current);
                }
                else
                {
                    np = null;
                }
            }
        }

        private void TrendlineBSpline(igpPointSet ptset)
        {
            Monomials.Clear();

            int count = ptset.Points.Count;

            if (count > 2)
            {
                IEnumerator item = ptset.Points.GetEnumerator();
                item.Reset();

                IgpPoint p1, p2, p3, p4;
                p1 = p2 = p3 = p4 = null;

                item.MoveNext();
                p1 = ((IgpPoint) item.Current);
                int it = 0;
                double xm = ((IgpPoint) item.Current).X;
                double tminx = xm;
                double tmaxx = xm;

                while (item.MoveNext())
                {
                    xm = ((IgpPoint) item.Current).X;

                    if (tminx > xm)
                    {
                        tminx = xm;
                    }

                    if (tmaxx < xm)
                    {
                        tmaxx = xm;
                    }

                    it++;

                    if (it == 1)
                    {
                        p2 = ((IgpPoint) item.Current);
                    }

                    if (it == ptset.Points.Count - 2)
                    {
                        p3 = ((IgpPoint) item.Current);
                    }

                    if (it == ptset.Points.Count - 1)
                    {
                        p4 = ((IgpPoint) item.Current);
                    }
                }

                DomainLower = tminx;
                DomainHigher = tmaxx;

                var main = new Matrix.Matrix(4*(count - 1), 4*(count - 1) + 1);

                //iteracja warunku zgodnosci wartosci //////////////////////////////////////////

                //przegladanie po punktach


                IgpPoint cp, np;

                int firstColumn = 0;
                int row = 0;

                item.Reset();

                item.MoveNext();

                cp = ((IgpPoint) item.Current);
                item.MoveNext();
                np = ((IgpPoint) item.Current);

                while (cp != null && np != null)
                {
                    main.Data[row, firstColumn] = Math.Pow(cp.X, 3);
                    main.Data[row, firstColumn + 1] = Math.Pow(cp.X, 2);
                    main.Data[row, firstColumn + 2] = cp.X;
                    main.Data[row, firstColumn + 3] = 1;

                    main.Data[row + 1, firstColumn] = Math.Pow(np.X, 3);
                    main.Data[row + 1, firstColumn + 1] = Math.Pow(np.X, 2);
                    main.Data[row + 1, firstColumn + 2] = np.X;
                    main.Data[row + 1, firstColumn + 3] = 1;

                    //wartosci w punktach
                    main.Data[row, 4*(count - 1)] = cp.Y;
                    main.Data[row + 1, 4*(count - 1)] = np.Y;

                    row += 2;
                    firstColumn += 4;

                    cp = ((IgpPoint) item.Current);
                    if (item.MoveNext())
                        np = ((IgpPoint) item.Current);
                    else np = null;
                }

                //iteracja warunku zgodnosci I i II pochodnej //////////////////////////////////

                //przegladanie po punktach

                item.Reset();

                item.MoveNext();
                item.MoveNext();
                cp = ((IgpPoint) item.Current);
                item.MoveNext();
                np = ((IgpPoint) item.Current);

                firstColumn = 0;

                while (np != null)
                {
                    //cp = ((IgpPoint)item.Current);
                    // I pochodna
                    main.Data[row, firstColumn] = 3*Math.Pow(cp.X, 2);
                    main.Data[row, firstColumn + 1] = 2*cp.X;
                    main.Data[row, firstColumn + 2] = 1;

                    main.Data[row, firstColumn + 4] = -3*Math.Pow(cp.X, 2);
                    main.Data[row, firstColumn + 5] = -2*cp.X;
                    main.Data[row, firstColumn + 6] = -1;

                    // II pochodna

                    main.Data[row + 1, firstColumn] = 6*cp.X;
                    main.Data[row + 1, firstColumn + 1] = 2;

                    main.Data[row + 1, firstColumn + 4] = -6*cp.X;
                    main.Data[row + 1, firstColumn + 5] = -2;

                    row += 2;
                    firstColumn += 4;

                    cp = ((IgpPoint) item.Current);

                    if (item.MoveNext())
                    {
                        np = ((IgpPoint) item.Current);
                    }
                    else
                    {
                        np = null;
                    }
                }

                // 2 rownania na spline naturalny

                main.Data[4*(count - 1) - 2, 0] = 3*Math.Pow(tminx, 2);
                main.Data[4*(count - 1) - 2, 1] = 2*tminx;
                main.Data[4*(count - 1) - 2, 2] = 1;

                if (p2.X/p1.X != 1)
                {
                    main.Data[4*(count - 1) - 2, main.Cols - 1] = (p2.Y - p1.Y)/(p2.X - p1.X);
                }

                main.Data[4*(count - 1) - 1, 4*(count - 1) - 4] = 3*Math.Pow(tmaxx, 2);
                main.Data[4*(count - 1) - 1, 4*(count - 1) - 3] = 2*tmaxx;
                main.Data[4*(count - 1) - 1, 4*(count - 1) - 2] = 1;

                if (p4.X/p3.X != 1)
                {
                    main.Data[4*(count - 1) - 1, main.Cols - 1] = (p4.Y - p3.Y)/(p4.X - p3.X);
                }

                Matrix.Matrix data = main.Gauss();

                int lastread = 0;
                item.Reset();
                item.MoveNext();

                cp = ((IgpPoint) item.Current);
                item.MoveNext();
                np = ((IgpPoint) item.Current);

                while (np != null)
                {
                    for (int y = 0; y < 4; y++)
                    {
                        IgpMonomial mon = AddMonomial(data.Data[lastread, 0], 3 - ((lastread)%4));
                        mon.LowerXLimit = Math.Min(cp.X, np.X);
                        mon.HigherXLimit = Math.Max(cp.X, np.X);
                        mon.Limited = true;
                        lastread++;
                    }

                    cp = ((IgpPoint) item.Current);

                    if (item.MoveNext())
                    {
                        np = ((IgpPoint) item.Current);
                    }
                    else
                    {
                        np = null;
                    }
                }
            }
        }

        /// <summary>
        /// Updates parameters of the trendline.
        /// </summary>
        public void TrendlineUpdate()
        {
            switch (TrendlineMode)
            {
                case IgpInterpolationMode.Linear:
                    TrendlineLinear(TrendlinePoints);
                    break;
                case IgpInterpolationMode.Lagrange:
                    TrendlineLagrange(TrendlinePoints);
                    break;
                case IgpInterpolationMode.BSpline:
                    TrendlineBSpline(TrendlinePoints);
                    break;
                case IgpInterpolationMode.BSplineNormal:
                    TrendlineBSplineNormal(TrendlinePoints);
                    break;
            }
        }
    }
}